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Reconstruction  of  processes  and  fields  from  noisy  data  is  to  solve  a  set  of  linear  algebraic  equa¬ 
tions.  Three  factors  affect  the  accuracy  of  reconstruction:  (a)  a  large  condition  number  of  the 
coefficient  matrix,  (b)  high  noise-to-signal  ratio  in  the  source  term,  and  (c)  no  a  priori  knowledge 
of  noise  statistics.  To  improve  reconstruction  accuracy,  the  set  of  linear  algebraic  equations  is 
transformed  into  a  new  set  with  minimum  condition  number  and  noise-to-signal  ratio  using  the 
rotation  matrix.  The  procedure  does  not  require  any  knowledge  of  low-order  statistics  of  noises. 
Several  examples  including  highly  distorted  Lorenz  attractor  illustrate  the  benefit  of  using  this 
procedure. 
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1.  Introduction 

Reconstruction  (reproduction  from  noisy  data)  of 
processes  and  fields  in  modern  physics,  geophysics, 
astrophysics,  plasma  physics  and  other  disciplinary 
sciences  is  to  solve  an  ill-posed  linear  algebraic 
equation, 


Aa  =  QY ,  (1) 

where  a  is  the  estimated  state  vector  (L- 
dimensional)  for  the  exact  state  vector  a;  A  is  a 
PxL  coefficient  matrix,  Q  is  a  P  x  P  square  matrix 
(P  >  L);  Y  is  a  P-dimensional  observation  vector, 
consisting  of  a  signal  Y  and  a  noise  Y', 

Y  =  Y  +  Y7. 

The  two  known  matrices  A  and  Q  are  determined 
by  the  physical  process  or  field.  Let  || . . .  ||  be  the 


Euclidean  norm  and 


IIQY'II 
IIQYII  ’ 


max(singular  values  of  A) 
min(singular  values  of  A) 

be  the  noise-to-signal  ratio,  dimension  ratio  and 
condition  number  of  the  matrix  A.  For  a  particular 
system,  772  is  given. 

Usually,  rji  and  773  are  large  (called  “imper¬ 
fect”), 

771  >  1,  773  >  1 , 

which  makes  (1)  difficult  to  solve.  Besides,  the  low- 
order  noise  statistics  and  the  norm  of  the  exact 
state  vector  (||a||)  are  unknown.  Reduction  of  771 
and  773  (i.e.  reduction  of  imperfection)  without 
a  priori  knowledge  of  noise  statistics  and  ||a||  is  an 
important  step  toward  solving  (1)  accurately.  If  the 


2991 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2QQ4  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2004  to  00-00-2004 

4.  TITLE  AND  SUBTITLE 

Rotation  Method  for  Reconstructing  Process  and  Field  From  Imperfect 
Data 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROTECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Ocean  Analysis  and  Prediction  Laboratory, Department  of 
Oceanography, Naval  Postgraduate  School, Monterey, CA, 93943 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

_ _ _  ABSTRACT 

18.  NUMBER  19a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  OS 

unclassified  unclassified  unclassified  Report  (SAR) 

7 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2992  P.  C.  Chu  et  al. 


matrix  A  has  no  noise,  the  accuracy  in  determining 
a  is  given  by  [Tikhonov  et  al.,  1990], 


2  _ 


all2 


<  vhi  ■ 


(2) 


Traditional  regularization  methods  (e.g. 
[Tikhonov  &  Arsenin,  1979;  Bennett,  1992])  trans¬ 
form  (1)  into  a  system  with  stable  solution  to  small 
perturbations  in  Q  and  Y.  However,  they  lead  to 
biased  estimations  a  for  a  and  do  not  guarantee 


a  — >  a  as  771  — >  0 . 

Furthermore,  the  traditional  regularization  meth¬ 
ods  require  a  priori  knowledge  of  ||a||  . 

To  overcome  these  weaknesses,  a  new  rotation 
method  for  772  <  1  is  developed  in  this  study  to 
change  (1)  into  a  new  system  with  possibly  min¬ 
imum  coefficient  matrix  and  noise-to-signal  ratio 
without  a  priori  knowledge  of  noise  statistics. 


2.  Rotation  Method 

When  the  dimension  ratio  772(<  1.0)  is  given,  large 
values  of  771  and  773  (called  “imperfect”  data), 

Vl  >  1)  V3  »  1 , 

make  (1)  difficult  to  solve.  Note  that  the  low-order 
noise  statistics  and  the  norm  of  the  exact  solution 
(||a||)  are  usually  unknown.  Reduction  of  771  and 
773  (i.e.  reduction  of  imperfection)  without  knowing 
noise  statistics  and  ||a||  is  an  important  step  toward 
solving  (1). 

Nonsingular  orthogonal  transformation  is  con¬ 
ducted  through  multiplication  of  (1)  by  a  plane 
rotation  matrix  S  from  the  left, 

SAa  =  SQY ,  (3) 


which  changes  the  coefficient  matrix  and  the  source 
term  from  (A,  QY)  to  (SA,  SQY)  and  provides 
the  opportunity  to  minimize  the  imperfection  of  the 
new  system  (3), 


where 


m 


viil+Vi) 


IISQY'H 
IISQYII  ’ 


— ►  min , 

(4) 

-  _  IISQYII 

INI  • 

(5) 

Minimization  (4)  is  obtained  by  the  following  max¬ 
imization  [Strakov,  1991;  Ivanov  et  al.,  2001] 


Ji  =  IIAf  - 


IISQYII1 


max, 


(6) 


where  ||A||  is  the  spherical  norm  of  the  matrix  A. 
Substitution  of  (5)  into  (6)  leads  to 


ji  =  IIaH2-»?3 


2(SQY*_SQYQ  2 
IISQYII* 


max, 

(7) 


which  is  the  procedure  to  obtain  minimum  values  of 
771  and  77 3  without  knowing  ||a]|2.  Here,  the  symbol 
indicates  the  scalar  product  in  the  Euclidean 
space.  For  a  white  noise  Y',  we  have 

(SQY*SQY')  0  as  P  — ►  00 . 

The  maximization  of  J\  using  (7)  is  then  equivalent 
to  the  minimization  (4).  Note  that  this  minimiza¬ 
tion  is  not  necessarily  the  same  as 

773  — >  min,  f/i  —►  min  . 


3.  Accuracy 


The  accuracy  of  the  reconstruction  (1)  is  estimated 
by  (2)  with  the  given  norm  of  the  exact  state  vec¬ 
tor  ||a||.  Since  ||a||  is  not  a  priori  known,  an  effective 
norm  ||aeff  ||  is  defined  by  the  following  minimization 
process, 


J2  = 


Ml* 

Null2 


— ►  min , 


(8) 


Simple  analysis  on  (4)-(8)  shows  that  two  norms 
Kffll  and  Kill  (generally  ||a*ff||  +  ||a*^|j)  exist 
that  the  functionals  J\  and  J2  reach  their  maximum 
values,  respectively,  and 


If  f|a||  - 
we  have 


INI  <  Kffll,  K5||<2||a||. 

Kffll  then  Hlf  0.  When  ||a|| 


*eff  I 


T  2 


1  - 


Ssffl 


Therefore,  replacement  of  ||a||  by  ||a*£||  will  not 
deteriorate  the  accuracy  of  the  solution  of  (1). 
A  similar  replacement  of  ||a||  is  also  suggested 
for  colored  noises  and  deterministic  perturbations. 
However,  the  reconstruction  accuracy  might  dete¬ 
riorate  because  the  scalar  product  (SQY  SQY') 
does  not  tend  to  0  for  P  — ►  00. 


4.  Applications 

The  new  transformed  system  (3)  can  be  solved  by 
usual  algebraic  methods  such  as  the  Gauss  method. 
For  simplicity,  we  restrict  the  reconstruction  to  data 
distorted  by  the  white  noise  only. 
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4.1.  Example  1 :  Linear  scalar 
process 

Consider  a  temporally  varying  linear  scalar  process, 

Y(t)=a  +  pt,  t  €  [0,  T] 

perturbed  by  a  white  noise  Y'(t).  The  noisy  data 
are  represented  by 

Y(t)  =  Y(t)  +  Y'(t)  . 


The  Fourier  series 


38 

X  ( t )  —  ao  + 

i=l 


+  at+19  cos 


(9) 


is  used  to  recover  the  signal  Y(t).  The  reconstruc¬ 
tion  is  to  determine  the  estimated  state  vector 


a  =  (ai,  a2, . . . ,  ai),  L  =  38, 

such  that  the  temporally  integrated  difference  be¬ 
tween  X  (t)  and  Y  (t)  reaches  the  minimum  value, 

rT 

J  =  /  [X(t)  -  Y(t)]2dt  -»  min  .  (10) 

Jo 

Substitution  of  (9)  into  (10)  leads  to  Eq.  (1)  of  a 
with  a  very  high  condition  number  for  the  coeffi¬ 
cient  matrix. 

The  noisy  dataset  Y(t)  is  generated  in  such  a 
way  that  the  sensitivity  of  reconstruction  quality  to 
(771 ,  772)  can  be  investigated:  60-100  points  (varying 
P  and  time  step  At)  uniformly  distributed  into  the 
interval  [0,  T]  are  used  to  get  Y(t),  and  white  noises 
Y'{t )  with  771  ranging  from  0  to  3  are  added  to  Y (t). 
Here,  T/ 60  >  At  >  T/100. 

The  three  nondimensional  parameters  771,  772 
and  773  vary  within  the  following  ranges, 

0  <  771  <  3;  0.4  <  772  <  0.6;  773  cs  106 .  (11) 

For  such  a  high  condition  number,  Eq.  (1)  cannot 
be  solved  directly.  The  rotation  matrix  S  is  calcu¬ 
lated  using  the  maximization  procedure  depicted  in 
Sec.  2.  The  condition  number  of  the  coefficient  ma¬ 
trix  of  the  new  system  (3)  reduces  to 

773  1.5 . 

The  reconstructed  accuracy  ||7||  reduces  with  the 
decrease  of  771  (noise-to-signal  ratio  before  the  rota¬ 
tion)  monotonically  from  4.0  for  771  =  3.0  to  0  for 


Fig.  1.  Reconstruction  of  the  one-dimensional  linear  pro¬ 
cess:  (a)  upper  panel:  reconstruction  accuracy,  and  (b)  lower 
panel:  comparison  between  Y  ( t )  (solid  line)  and  recon¬ 
structed  processes  X(t)  (dashed  curves)  with  1,  2,  3,  4  and  5 
corresponding  to  771  =  0.6,  1.2,  1.8,  2.4,  and  3.0. 


771  =  0  [Fig.  1(a)].  The  reconstructed  process  X(t) 
is  closer  to  Y(t)  as  771  reduces  [Fig.  1(b)]. 

4.2.  Example  2:  Two  dimensional 
field 

Eremeev  et  al.  [1992]  reconstructed  the  Black  Sea 
summer  climatological  surface  temperature  field 
(TcUm)  using  the  traditional  regularization  method 
[Bennett,  1992].  The  nondimensional  parameters 
(771,  772,  773)  vary  within  the  following  intervals: 

0  <  771  <  4,  0.025  <  772  <  0.67, 

3  x  104  <  773  <  3  x  107 . 

For  large  values  of  (771,  773),  the  accuracy  of  the  tra¬ 
ditional  regularization  methods  is  not  very  good. 

Consider  that  the  climatological  surface  tem¬ 
perature  field  Tciim(:r,  y)  is  perturbed  by  white 
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Fig.  2.  Dependence  of  the  accuracy  ||7||  on  771  with  different  revalues  for  the  two-dimensional  temperature  field  in  the  Black 
Sea.  Note  that  ||-y||  is  usually  smaller  than  y  1. 


noises  T'(x,  y).  The  noisy  data  are  represented  by 
TdataC^j  V )  =  Tciim(x,  y)  +  T  ( x ,  y),  {x,  y)  €  ft  , 

where  ft  is  the  Black  Sea.  The  generalized  Fourier 
series, 

30 

T(x,y)  =  Y^aiMx,  V ) .  (13) 

i=l 

is  used  to  recover  the  signal  Tciim(x,  y).  Here,  the 
basis  functions  {^i(x ,  y),  i  =  1,  2, ... ,  30}  are  the 
eigenfuctions  of  the  plane  von  Neumann  operator 
with  homogeneous  boundary  conditions  for  the  do¬ 
main  ft  [Eremeev  et  al.,  1992].  The  reconstruction 
is  to  determine  the  estimated  state  vector 

a  =  (oi,  a2, . . . ,  ai),  L  =  30, 

such  that  the  spatially  integrated  difference  between 
Tdata  (x,  y)  and  T(x,  y)  reaches  the  minimum  value, 

J  =  [[  [T(x,  y)  -  Tdata (x,  y)]2dxdy  -*  min  .  (14) 
J  Jn 

Substitution  of  (13)  into  (14)  leads  to  Eq.  (1)  of  a 
with  a  very  high  condition  number  for  the  coeffi¬ 
cient  matrix  [see  (12)]. 

After  applying  the  rotation  matrix  S  with  the 
maximization  procedure  [Eq.  (7)],  the  condition 


number  of  the  coefficient  matrix  of  the  new  system 
(3)  reduces  to 


773  <  5.0 . 

The  reconstructed  accuracy  ||7||  of  our  recon¬ 
struction  scheme  is  up  bounded  by  771  (noise-to- 
signal  ratio  before  the  rotation) 

Ibil  <  Vi , 

and  usually  decreases  with  decreasing  771  and  772 
(Fig.  2).  Note  that  the  approach  does  not  improve 
the  reconstruction  accuracy  for  772  <  0.045.  For 
large  noise-to-signal  ratio  (771  >  3.0),  the  recon¬ 
struction  accuracy  improves  when  the  dimension 
ratio  772  increases  from  0.3  to  0.67. 


4.3.  Example  3:  Lorenz  Attractor 

Lorenz  system,  a  truncated  three-component  atmo¬ 
spheric  convection  model,  is  represented  by  a  three- 
dimensional  vector 
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which  satisfies  the  following  nonlinear  ordinary 
differential  equation  with  the  initial  condition 

^  =  F(X,  a),  X(to)  =  X0  ,  (15) 

at 

where  the  vector  functional  F  and  state  parameter 
vector  are  defined  by  [Lorenz,  1963] 


— aiXi  +  0,3X2 

F  = 

—X1X2  +  02X1  —  X2 

,  a  — 

02 

X1X2  —  03X3 

-a3. 

The  initial  condition  (X0)  and  the  parameter  vector 
(a)  affect  the  characteristics  of  the  Lorenz  system. 
Chu  [1999]  showed  that  boundary  condition  was 
represented  by  the  model  parameter.  We  integrate 
Eq.  (15)  with  the  time  step  At  =  0.01  and 


'10' 

1 - 

7— t  7-H 

0  d 
■ 

,  a  ~ 

28 

L0.1. 

8 

.  3  . 

to  obtain  the  famous  butterfly  pattern  for  the  track 
in  the  phase  space  (Xi,  X2)  [Fig.  3(a)]  and  unstable 
oscillation  for  temporal  variation  of  X3  [Fig.  3(b)]. 

The  observation  vector  of  the  Lorenz  sys¬ 
tem  [Y(£)]  contains  signal  [Y(t)}  and  white  noises 

[Y  '(£)], 

Y(£)=Y(£)  +  Y'(t) 

with  7 h  =  0.6.  The  butterfly  pattern  is  totally  de¬ 
stroyed  in  the  phase  space  (Yi,  >2)  [Fig.  4(a)].  The 
time  series  of  the  third  component  (Y3)  shows  a 
stochastic  process^ [Fig.  4(b)]. 

The  signal  Y(f)  is  recovered  from  the  noisy 
data  [Yi(f),  Y2(t),  Y3(£)]  (Fig.  4)  through  accurate 
estimation  of  the  parameter  vector  a  and  effective 
reduction  of  noises.  Estimation  of  the  state  vector 

a=  (ai,. ,  <xl),  L  =  3, 

is  conducted  using  the  minimization  of  the  tempo¬ 
rally  integrated  difference  between  X(£)  and  Y (£), 

J  =  [  fX(t)  -  Y(t)]*[X(«)  -  Y (t)]dt  -»  min  . 

JtQ 

(18) 


Different  from  Example  2,  there  is  no  explicit  re¬ 
lationship  between  X(t)  and  a.  The  minimization 
(18)  becomes 


dJ 

daj 


=  0. 


An  iterative  algorithm  [Eykhoff,  1973]  and  the 
initial  condition 

X(t0)  =  Y(£o)  ,  (20) 


are  used.  Let  the  parameter  vector  be  estimated  by 
a(n_1)  for  the  ( n  -  l)th  iteration.  Integration  of  the 
Lorenz  system 


0.0  — ^ - 1 - 1  |  I  \* 

0.0  5.0  10.0  15.0  20.0  25.0  30.0 

(b) 

Fig.  3.  Lorenz  attractor  without  white  noise:  (a)  butter¬ 
fly  pattern  of  track  in  (Xi,  X2)  phase  space,  (b)  time  series 
of  A3. 


(19) 
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0.0  5.0  10.0  15.0  20.0  25.0  30.0 


(b) 

Fig.  4.  Lorenz  attractor  distorted  by  white  noise  with  rji  = 
0.6:  (a)  noisy  track  in  (Yi,  Y2)  phase  space  with  no  butterfly 
pattern,  (b)  time  series  of  Y3. 


0.0  5.0  10.0  15.0  20.0  25.0  30.0 


(b) 

Fig.  5.  Reconstructed  Lorenz  attractor  from  the  “noisy 
data”  with  771  =  0.6  using  the  rotation  method:  (a)  butter¬ 
fly  pattern  of  track  in  (Yi,  Y2)  phase  space,  (b)  time  series 
of  Y3. 


gives  the  state  vector  for  the  (n  —  l)th 

iteration.  If  the  estimated  parameter  vector  a/"-1) 
has  an  increment  Aa,  the  state  vector  has  a  corre¬ 
sponding  increment, 

AX(n— i)(t)  =  u(n_1)Aa,  (22) 

where  U  is  the  sensibility  matrix  defined  by  Uij  = 
dXi/daj;  i,  j  =  1,  2,  3.  When  the  increment  Aa 


is  selected  such  that  Eq.  (19)  is  satisfied,  this  leads 
to 

T 

f  (U*)!"-1)!^"-1*  *  Aadt 

Jto 

=  [T *  ( Y  —  X(n_1))dt .  (23) 

Jto 

The  temporally  integrated  difference  between  Y  ( t ) 
and  X(n_1)(t)  +  TJ(n-1)Aa,  reaches  the  minimum 
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[see  (18)].  Here,  the  matrix  U4  is  the  transpose  of 
the  matrix  U.  The  estimated  parameter  vector  at 
the  nth  iteration  is  taken  as 

a .(")  =  a(n_1)  +  Aa .  (24) 

We  solve  Eq.  (21)  with  a*”)  and  get  the  state  vector 
at  the  nth  iteration,  X(n)(t). 

At  each  iteration  step,  an  ill-posed  algebraic 
equation  for  Aa  [Eq.  (23)]  should  be  solved.  As  for 
the  two  previous  examples,  the  noise-to-signal  ratio 
and  the  condition  number  axe  large, 

7/1  ~  0.6,  7/3  ~  104. 

The  rotation  method  is  used  to  get  minimum  (7/1, 
7/3)  at  each  iteration  step.  The  ill-conditioning  of 
Eq.  (23)  is  greatly  reduced  because 

m  <  5, 

at  all  iteration  steps. 

The  iteration  converges  fast.  The  final  esti¬ 
mated  parameter  vector  is  given  by 

9.93  \ 

27.98  j  . 

2.667/ 

After  the  reconstruction,  we  obtain  the  butterfly 
pattern  of  the  track  in  the  phase  space  (Xi,  X2) 
[Fig.  5(a)]  and  the  unstable  oscillation  for  tempo¬ 
ral  variation  of  X3  [Fig.  5(b)].  Comparison  between 
Figs.  3  and  5  demonstrates  the  capability  of  the 
rotation  method  for  reconstructing  noisy  data. 

5.  Conclusions 

(1)  The  rotation  method  has  two  major  advantages 
in  reconstruction  of  process  and  field  from  imperfect 
data.  First,  it  does  not  need  any  a  priori  knowledge 
of  noise  statistics  and/or  representation  for  deter¬ 
ministic  perturbations.  Second,  it  reduces  the  con¬ 
dition  number  and  noise-to-signal  ratio  drastically, 
and  in  turn  reduces  the  imperfection  of  the  system. 

(2)  The  rotation  method  has  capability  to  effec¬ 
tively  reconstruct  the  linear  processes  with  ranges 
of  0  <  7/1  <  4  —  5;  0  <  772  <  0.5  —  0.6;  7/3  <  106, 
and  the  nonlinear  Lorenz  attractor  with  ranges  of 
0  <  Vi  <  0-6;  0  <  7/2  <  0.5  —  0.6;  7/3  <  104  at  each 
iteration. 

(3)  The  procedure  can  be  applied  to  processes 
and  fields  with  both  stochastic  and  determinis¬ 
tic  perturbations.  Behavior  of  the  scalar  product 


(SQY*SQY')  —7  0  determines  the  effectiveness  of 
noise  reduction.  Besides,  it  can  be  easily  general¬ 
ized  on  stochastic  signals  if  there  is  no  correlation 
between  Y  and  Y' . 
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